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Abstract 

This paper introduces a new four-parameter lifetime model called the Weibull 
Birnbaum-Saunders distribution. This new distribution represents a more 
flexible model for the lifetime data. Its failure rate function can be increas¬ 
ing, decreasing, upside-down bathtub shaped, bathtub-shaped or modified 
bathtub shaped depending on its parameters. Some structural properties of 
the proposed model are investigated including expansions for the cumulative 
and density functions, moments, generating function, mean deviations, or¬ 
der statistics and reliability. The maximum likelihood estimation method is 
used to estimate the model parameters and the observed information matrix 
is determined. The flexibility of the new model is shown by means of two 
real data sets. 


Keywords: Birnbaum-Saunders distribution, Weibull-G class, moment, or¬ 
der statistic, maximum likelihood estimation, observed information matrix. 


1 Introduction 

The two-parameter Birnbaum-Saunders (BS) distribution, which was intro¬ 
duced by Birnbaum and Saunders [3], is a very important lifetime model and 
is widely used in reliability studies. This distribution, also known as the 
fatigue life distribution, was originally derived from a model that shows the 
total time that passes until that some type of cumulative damage, produced 
by the development and growth of a dominant crack, surpasses a threshold 
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value and causes the material specimen to fail. Desmond [10] provided a 
more general derivation based on a biological model and strengthened the 
physical justification for the use of this distribution. 

A random variable T having the BS distribution with shape parameter a > 0 
and scale parameter (3 > 0, denoted by T ~ BS(a,/3), is defined by 


T = (3 
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where Z is a standard normal random variable. The cumulative distribution 
function (cdf) of T is 

G(t) = $(v), t> 0, (1) 

where $ (•) is the standard normal distribution function, v = a~ 1 p ( t//3 ) and 
p (z) = z 1//2 — z~ 1 ^ 2 . The probability density function (pdf) corresponding to 
(HD is given by 

g(t) = K(a,p)t~ 3/2 (t + /3)exp j , t > 0, (2) 


where k (ot,/3) = exp (a 2 ) / (2a\/27r/3) and r (z) = z + z h The fractional 
moments of T are given by (see |20|) 

E(T i )=P t I(a,P), 


where 


I ( a,/3 ) 


Kk+ 1/2 (« 2 ) + Kk-l/2 (« 2 ) 

2Ad/ 2 («- 2 ) 


(3) 


and the function K v ( z ) denotes the modified Bessel function of the third 
kind with v representing its order and z the argument. The parameter (3 is 
the median of the BS distribution, because G (f3) = (0) = 1/2. 


Since the BS distribution was proposed, it has received much attention in 
the literature. This attention for the BS distribution is due to its many 
attractive properties and its relation to the normal distribution. For more 
details on the BS distribution, we refer to 15] and the references therein. The 


BS distribution has been used in several research areas such as engineering, 
environmental sciences, finance, and wind energy. However, it allows for 
upside-down hazard rates only (see 0), hence cannot provide reasonable 
fits for modeling phenomenon with decreasing, increasing, modified bathtub 
shaped and bathtub-shaped failure rates which are common in reliability 
studies. 
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For this reason, several generalizations and extensions of the BS distribution 
have been proposed in the literature. For example, Cordeiro and Lemonte 
[6], using the beta-G class fill ] . proposed an extension of BS distribution 
named as the beta BS distribution. Sanlo et al. [22], based on the work 
of Cordeiro and de Castro [iij], defined the Kumaraswamy BS distribution. 
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Lemonte 16], based on the scheme introduced by Marshall and Olkin 
deihned the Marshall-Olkin extended BS distribution. Cordeiro et al. 
adopted the McDonald-G class [2] to dehne the McDonald BS distribution 
Cordeiro et al. [2] used the generator approach of Zografos and Balakrishnan 
[28] to introduce the gamma Birnbaum-Saunders distribution. In this paper, 
a new four-parameter extension for the BS distribution is proposed. 


Recently, Bourguignon et al. [4] proposed an interesting method of adding 
a new parameter to an existing G distribution. The resulting distribution, 
known as the Weibull-G distribution, gives more flexibility to model various 
types of data. Let G (t, 9) be a continuous baseline distribution with density g 
depends on a parameter vector 9 and the Weibull cdf F w (w) = 1 — e " (for 
w > 0) with positive parameters a and b. Bourguignon et al. [ij replaced the 
argument w by G (w,£) jG ( w , 9) where G (u>, £) — 1 — G (w, 9), and defined 
the cdf of their class by 


r c(M) 

F (f; a,b,9) = ab j ( ’ )J w b ~ 1 e~ awh dw = 1 — exp ( —a 


G (f, 9) 
_G (t, 9)_ 


( 4 ) 


Then, the Weibull-G density function is given by 

'G(t,9) b ~ 1 ' 


f (t; a, b, 9) = abg (t, 9) 


G (t, 9) b+1 


exp —a 


G(t,9) 
G (t, 9) 


( 5 ) 


Some generalized distributions have been proposed under this methodology. 


Tahir et al. [23|, [24J, l25j] dehned the Weibull-Pareto, Weibull-Lomax and 


Weibull-Dagum distributions by taking G (t , 9) to be the cdf of the Pareto, 
Lomax and Dagum distributions, respectively. More recently, Ahfy et al. [lj 
dehned and studied the Weibull Frechet distribution. In a similar way, we 
propose a new extension for the BS distribution called the Weibull BS (WBS) 
distribution, which has been applied to the modeling of fatigue failure times 
and reliability studies. 


The rest of the paper is organized as follows. In Section [2] we introduce the 
WBS distribution and plot its density and failure rate functions. In Section 
[31 we provide a mixture representation for its density and cumulative dis¬ 
tributions. Structural properties such as the ordinary moments, generating 
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function, quantile function and simulation, mean deviations, the density of 
the order statistics and the reliability are derived in Section U In Section 
[5l we discuss maximum likelihood estimation of the WBS parameters and 
derive the observed information matrix. Two applications are presented in 
Section to show the potentiality of the new distribution. Some concluding 
remarks are given in Section [3 


2 The WBS distribution 


Substituting (HD in (HD, the cdf of the WBS distribution can be written as 


F(t) 


1 — exp 



*(v) ] b \ 

l-$(u)J J 


( 6 ) 


The pdf corresponding to (ED is given by 

/ ( t) = abn (a, (3) t~ 3 ^ (t + /3) exp 
< f > 0) 6 " 1 


(*/m 


X 


{!-*(«)} 


6+1 


exp —a 


2a 2 J 

$(t>) 

1 - $ (v) 


(7) 


where j3 is a scale parameter and a, a and b are positive shape parameters. It 
is clear that the BS distribution is not a special case of WBS distribution. If 
a random variable T follows a WBS distribution with parameters a, (3, a and 
b will be denoted by T ~ WBS(ct, f3, a, b). The reliability and the failure rate 
function of T are, respectively, given by 


R{t) = exp 



!-$(n)J J 


) 


and 


h(t) = abn ( a , (3) t 3 ^ 2 (t + (3) 


(v) 

{!-»(«))“ 


exp 


2 a 2 /' 


Plots of pdf and failure rate function of the WBS distribution for selected 
values of the parameters are given in Figure 1 and Figure 2, respectively. Fig¬ 
ure 1 indicates that the WBS pdf can take various shapes such as symmetric, 
right-skewed and left-skewed depending on the parameter values. Figure 2 
shows that the failure rate function of the WBS distribution can be increas¬ 
ing, decreasing, upside-down bathtub (unimodal) shaped, bathtub-shaped or 
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modified bathtub shaped (unimodal shape followed by increasing) depending 
on the parameter values. So, the WBS distribution is quite flexible and can 
be used effectively in analyzing survival data. 




t t 

Figure 1. Plots of the WBS pdf for some values of the parameters. 




t t 

Figure 2. Plots of the WBS failure rate function for some values of the 

parameters. 
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3 Mixture representation 


In this section, we derive expansions for the pdf and cdf of the WBS dis¬ 
tribution. The pdf and cdf of the WBS distribution can be written as a 
linear combination of the pdf and cdf of exponentiated BS (EBS) distribu¬ 
tion, respectively. A random variable X having the EBS distribution with 
parameters a, (3 and a > 0, denoted by X ~ EBS(a, (3, a), if its cdf and pdf 
are given by H ( x ) = (u) and h(x ) = ag(:r) < I >a ~ 1 (T), respectively, where v 

is defined in (JT]) and g is given in ([2]) . The properties of exponentiated distri¬ 
butions have been studied by several authors. For example, Mudholkar and 
Srivastava [l9| studied the exponentiated Weibull distribution, Gupta and 
Kundu [l3[ studied the exponentiated exponential distribution and Sarhan 
and Apaloo [21] proposed the exponentiated modified Weibull extension dis¬ 
tribution. 


The pdf of the WBS distribution (JT]) can be written as 


f (t) = abg ( t ) 


$ 


6-1 


a) 


16+1 


exp —a 


$(v) 

1 — $ (v) 


[1-*W 

Using the series expansion for the exponential function, we obtain 

b\ °° / i\fc fc Abfc 


( 8 ) 


exp —a 


Substituting 

f (t) = abg (t) Y 



^(-1 )a k <h bfe (u) 

^ k\ [1 - $(n)] 6fc ' 


0) [1 - $ (u)n (6fc+?,+1) 


(9) 


k =0 


Since 0 < (v) < 1, for t > 0 and ( bk + 6 + 1) > 0, then by using the 

binomial series expansion [1 — $ (n)] _ bfe+b+1 given by 

(6fc+6+i) _ T (bk + b + 1 + j) 


[1 - $ (r,)] 

where T (•) is the complete gamma function, we obtain 


+- 1 ^ 1 ^ 1 ( v ) 

j^ o j\T(bk + b + 1) K h 


where 


/ (k) 'y Wkjhbk+b+j (k ), 

k,j =o 

(—l) fe ba k+1 T (bk + b+ 1 + j) 

Wk ’ j ~ k\j\(bk + b +j)Y(bk + b + iy 


( 10 ) 
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and h bk+b+ j (t) denotes the EBS(cq (3,bk + b + j ) density function. By inte¬ 
grating (fTUJ) . we get 


F{t)= ^ w kJ $ bk+b+j (v ). (11) 

k,j=0 

It easy clear that Ylh=o^%o Wk d = E Equation (TTOh means that the pdf 
of the WBS distribution is a double linear mixture of the pdf of EBS dis¬ 
tribution. Based on this equation, several structural properties of the WBS 
distribution can be obtained by knowing those of the EBS distribution. For 
example, the ordinary, inverse and factorial moments, generating function 
and characteristic function of the WBS distribution can be obtained directly 
from the EBS distribution. 

If b is a positive real non-integer, we can expand (x;) as 

OO 

Qbk+b+j ^ _ 

(bk + b + j)$ r (v), (12) 

r =0 

where 



Thus, from equations ([2]), (HUD and (HUD . we get 

OO 

f(t) = 9 (*) dr ® r ( v ) > ( 13 ) 

r =0 


where 

OO 

d r = u>h,jS r (bk + b + j). 

k,j =o 

4 Some structural properties 

In this section, we give some mathematical properties of the WBS distribu¬ 
tion. 

4.1 Moments 

In this subsection, we derive the expression for sth order moment of WBS 
distribution. The moments of some orders will help in determining the ex¬ 
pected life time of adevice and also the dispersion, skewness and kurtosis in a 
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given set of observations arising in reliability applications. The sth moment 
of the WBS random variable T can be derived from the probability weighted 
moments of the BS distribution. The probability weighted moments of the 
BS distribution are formally defined, for p and r non-negative integers, by 

poo 

Tp,r = / t p g (t) <f> r (v) dt. (14) 

Jo 

There are many softwares such as MAPLE, MATLAB and R that can be used 
to compute the integral (TT4|) numerically. From Jfjj], we have an alternative 
representation to compute T p<r that is 



where Sj — k\ + ... + kj, A (ki ,..., kj ) = a 2sj 3 (ik, ... akj, 
a*, = (—l) fc 2 ( J~ 2k )/ 2 (2 k + 1) k\] 1 and I (p + (2 Sj + j — 2m) /2, a) is 

calculated from ((2]) in terms of the modified Bessel function of the third 
kind. 


Therefore, the sth moment of T can be written from (TTH1) as 

OO 

d's = ^ d rT s ,r, (16) 

r =0 

where t S)T is obtained from (TT5|) and d r is given by (TT3|) . We can compute 
numerically the sth moment in any symbolic software by taking in the sum 
a large number of summands instead of infinity. 


4.2 Moment generating function 


Let T ~ WBS(a, /3, a, b ). The moment generating function of T, say M (s) = 
E (e sT ), is an alternative specification of its probability distribution. Here, 
we provide two representations for M (s). From ((131) . we obtain a first rep¬ 
resentation 

OO /*C>0 

M (s) = V d r / e st g (t) (u) dt. 

r=0 

Expanding the exponential function, we can rewritte M (s) as 


m w = EE 

r —0 p =0 




PI 


p ’ r s p . 




The second representation forM (s) is based on the quantile expansion of the 
BS distribution. From fflU]) . we have 



e st g (t) $ bk+b +J- 1 (v ) dt, 


where g(t) is the BS(a,/3) pdf. By setting u — T (v) in above integral, we 
get 


OO -l 

M (s) = Wk,j / u bk+b+:i ~ 1 exp (sQ (u)) du, 

k,j=o 


where t = Q(u ) is the quantile function of the BS distribution and u — $ (v) 
is given by (PQ). Using the exponential expansion, we get 


OO OO p 

M (s) = ]T w Ki h « 6fc+b+i- V («) du. (17) 

k,j =0 p=0 P' do 

From [&], if the condition —2 < (t/fd) 1 ^ 2 — (t//d)~^ 2 < 2 holds, we have the 
expansion for the quantile function of the BS distribution 

OO 

t = Q («) = ^2 Vi {u - 1/2)*, (18) 

i=0 

where 

OO 

Vi = (2vr)* /2 ^2 

3=0 

do = (3, d 2 j+i = / 5a 2 - ?+1 (V 2 )2” 2j for j >0,d 2 = /3a 2 /2, d 2 j = 0 for j > 2 and 
the quantities e 3/l can be determined from the recurrence equation 


i 

6j i (ia 0 ) ^ ^ (mj ~t~ m i) CL m Cj^— m , 

m =1 


and ej t o = a J 0 . Here, the quantities a m are defined by a m = 0 (for m = 
0, 2,4,...) and a m = 6( m _i)/2 (for i = 1, 3, 5,...), where the b m 's are com¬ 
puted recursively from 


bm+l — 


---V 

2 (2m+ 3) ^ 

\ / r= i 


(2 r + 1) (2m — 2 r + 1) b r b m _ r 
(r + 1) (2r + 1) 


The first constants are We have b 0 — 1, b\ — 1/6, b 2 = 7/120, bo = 
127/7560,.... 
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Inserting (|IS|) in f[T7jl . we get 


M ( s ) = 


OO OO p pi / oo 

E”^Ey/““ E 

fe,j=0 p=0 1/0 \ i=0 


rji (u — 1/2)* ] du. (19) 


From 12, Sec. 0.314] for a power series raised to a positive integer p, we have 


5>(«-i/2)* = 5^ fa - v 2 ) 1 > 


p=o 


t=0 


where the coefficients 5 Pi j (for i = 1,2,...) can be determined from the recur¬ 
rence equation 


8 P ,i = (ffio) 1 ^2 ( m P + m ~ i ) 


m=l 


and 5 Pt o = a? Q . Hence, comes directly from <5 Ji0 ,..., Sj t j_i and, therefore, 
from a 0 ,..., a*. Then 


00 ^ S P r 1 

M ( s ) = "^2 w k,j—^^p,i 

k,j,p,i= 0 ^ 


W 


6(fc+l)+j-l ( u _ X / 2 )i du. (20) 


Therefore, using the binomial expansion in (|2U1) . we obtain 

- 1 y~ l { i i) w k,A 


00 i 


m(°)= E E 


™— s r. 


k,j,i,p =0 1=0 


p\ (bk + b + j + l) 2 l ~ l 


4.3 Quantile function and simulation 

In this subsection, we give an expression for WBS quantile function, Q (u) = 
F _1 (rx), in terms of the BS quantile function Q B $ (•) . The BS quantile func¬ 
tion is straightforward computed from the standard normal quantile function 
<h _1 (u). We have (see , |.6]) 


Qbs fa) = J fad* 1 fa) + \A + [a®- 1 fa)] 2 ^ . 
Then, by inverting F(x) = u, we obtain 

Qwbs fa) = | fa) + ^4+ [a^- 1 fa)] 2 ^j , 
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where 


1+ (-i ln ( 1_M ))' 

Therefore, it is easy to simulate the WBS distribution. Let U be a contin¬ 
uous uniform variable on the unit interval (0,1]. Thus, using the inverse 
transformation method, the random variable T given by 

T = Qwbs{U ) = | (W 1 (P) + ^4 + (P)] 2 ) , (21) 


where 

p (->(1 

i+(->(i-t/)h’ 

has the WBS distribution. Equation (12TT) may be used to generate random 
numbers from the WBS distribution when the parameters are known. 

We plot the exact and the empirical cdf of WBS distribution in Figure 3 using 
a pseudo random sample of sizelOOO to check the correctness of the procedure 
for simulating a data set from WBS distribution. The histograms for two 
generated data sets and the exact WBS density plots from two simulated data 
sets for some parameter values are given in Figure 4. These plots indicate 
that the simulated values are consistent with the WBS distribution. 



Figure 3. Comparison of exact and empirical cdf of the WBS 
distribution to simulate random numbers. 
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(b) 


(a) 




0.5 1.0 1.5 2.0 2.5 3.0 3.5 


0.2 0.4 0.6 0.8 


Figure 4. Plots of the WBS densities for simulated data sets: (a) 
a = 2.5, /3 = 2.5, a = 2 and 6 = 4; (b) a = 0.5, (3 = 1.7, a = 0.2 and 

6 = 0.4. 


4.4 Mean deviations 

Let T be a random variable having the WBS(a, /?, a, 6) distribution. The 
mean deviations of T about the mean and about the median can be used as 
measures of spread in a population. They are given by 

poo 

5 1 = E{\T-^\)= \T — //jJ f (t) dt 
Jo 

and 

poo 

62 = E (|T — m\) = / \t — m\ f (t) dt, 

Jo 

respectively, where the mean fi\ is calculated from (1161) and the median m 
is given by m = Q WBS (l/2). The measures hi and h 2 can be expressed as 

hi = 2 n\F (//J — J (ni) and h 2 = E (|T — m\) = fi\ — 2 J (m) , 

where J (q ) = j^ tf (t) dt. From ([15]) . J ( q) can be written as 

OO 

J (?) = dri p r ) > ( 22 ) 

r =0 
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where 

<P (q, r )= I *9 (t) (v) dt - 

Jo 

From [g|, we have 

k (a, /3) ( r N 


ip(q, r ) = 


2 Sj -\-j 


E h V p-v.,+m A (tl ,..., kj) V (-py 


j =1 x ' k\,...,kj=0 


X 


( 2S ^ ^ J\(2s j+j -2m-l)/2 {t + ^ 


exp 


m=0 

T tt/P) 


2a 2 


dt. 


where Sj and A (ky ,..., kj) are defined in fTTHl) . Consider 


D (p, q) = / exp 


(t//3 + p/t) 


2 a 2 


rq/P 


dt = (3 p+1 / u p exp 


n + w 


-l 


2a 2 


du. 


From 26], we can write 


D(p,q) = 2P p+1 K p+ i (a~ 2 ) - g p+I K p+l ( 


q P 


2a 2 /?’ 2ct 2 g/ 


where, K v ( 21 , 22 ) is the incomplete Bessel function with order v and argu¬ 
ments 21 and 22 . Then, we obtain 


<p(q,r) = 


k (a , j3) r 


E 

i=i 


x 


2 s 7 + j 


m 


J 

D 


2 Sj + j 


E r p ‘‘ +i)/2 A(k 1 ,...,k j )'£(-P)’ 

ki,...,kj=0 

2Sj + 3 2m + l d 


m= 0 


2 s,- + j — 2m — 1 


which can be calculated from the function D ( p,q ). Hence, we can use this 
expression for <p(q,r) to compute J (q). From (l22lh we obtain the Bonfer- 
roni and Lorenz curves defined by B (p) = J (q) /p/i.\ and L(t) — J ( q ) //ij, 
respectively. These curves have applications in reliability. 


4.5 Reliability 

In the stress-strength modelling, R = P(T 2 < Tj) is a measure of compo¬ 
nent reliability when it is subjected to random stress To and has strength 
Tj. The component fails at the instant that the stress applied to it ex¬ 
ceeds the strength and the component will function satisfactorily whenever 
T\ > T 2 . The parameter R is referred to as the reliability parameter. This 
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type of functional can be of practical importance in many applications. In 
this Section, we derive the reliability R when T\ and T 2 have independent 
WBS(ct, /3, a\, b \) and WBS(cp /3, a 2 , b- 2 ) distributions. The pdf of 7\ and the 
cdf of T 2 can be obtained from (TTOT) and (TTTT) as 

OO 

k,j=0 


and 


respectively, where 


F 2 ( t ) = ^ w 2l , m ^ l+h 

1 , 171 =0 


'2 +m 


and 


We have 


Then 


where 


w lk,j ~ 


W21,m 


{-l) k b l a'i +1 T{b 1 k + b l + 1+j) 

k\j\r (61 k + 61 + 1 ) 

(—1) 6 2 a2 +1 T (b 2 l + b 2 + 1 + m) 
l\m\ ( b 2 l + b 2 + m)T ( b 2 l + b 2 + 1)' 

poo 

R= / fi(t)F 2 (t)dt. 


R= ^2 w* lk j w 2hm / g (t) <l> d * (■ V ) dt , 

k,j,l,m =0 *'° 


d* — fei (fc + 1) + b 2 (7 + 1) + j + 777 — 1. 
From (fT2]i , we can write 


* <! » = I>rW* r (»), 

r=0 


and then, we get 

OO OO 

d? ^ ^ 1 k.j ^' 2 /. m ^ ^ r (h) To^) 

k,j,l,m =0 r=0 

where r 0ir is obtained from (TT5]i . 
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4.6 Order statistics 


In this section, the distribution of the itii order statistic for the WBS distri¬ 
bution are presented. The order statistics play an important role in reliability 
and life testing. Let Ti,..., T n be a simple random sample from WBS distri¬ 
bution with cdf and pdf as in (j7D and (ED, respectively. Let T l n < ■ ■ ■ < T n n 
denote the order statistics obtained from this sample. In reliability literature, 
the ftli order statistic, say T i;n , denotes the lifetime of an (n — i + l)-out-of-n 
system which consists of n independent and identically components. The 
pdf of Ti- n is given by 


fi,n (t) = 


n! 




[n 


1 )! ^ 

' 3=0 


‘I ) J 


n — i 
3 


f (t) F it) 


i+j —1 


(23) 


From (ED, we have 

F i+j ~ l (t) = 

Using the binomial series expansion, we get 

F ,:+J - 1 (t) = jr (-1)"' ( l+ {~ exp ( -ka 


( 

\ (v) ] 

b )l 

1 — exp —a 

1 — $ (v) 

J 


i+j -1 


k =o 


k 


$(u) 

1 - <3> (v) 


(24) 


Inserting (JHD and (124p in (12U1) . we obtain 
n\abn (a, f3) f _3/2 (t + /3) 

fi,n (t) = - — - 77777 - 77 -exp 


(n — i)\ (n — 1)! 

n—k oo 


f T(t//9)) 

$(u) b 1 

l 2cU J 



x EE(-d 

j =0 k =0 


k+j 


i+j — l\fn —i 


k 


3 


exp —a (k + 1) 


$(u) 

1 - <3> (u) 




Using the power series for the exponential function, we have 


exp —a (k + 1) 


$(t>) 

1 — $ (v) 


(—l) z (ak + a) 1 $ ( v) bl 

h 1 [i - «(«-)]‘ 


Then 


fi,n (t) = 


g(t)n!ab 


n—k oo 


EEH 


k+j 


(n — i)\ (n — 1)\ / — 
y ’ v ’ j =o k =o 


i + j — 1\ (n — i 


x 


J2 (vf+h- 1 [1 - $ (v)]~ {bl+b+1) . 


1=0 
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Since 0 < $ (v) < 1, for t > 0 and (bl + b + 1) > 0, we have 




m =0 


Therefore, the pdf of the ith order statistic for WBS distribution is 


fi,n (t) ^ ^ $l,mhbl+b-\-m (t') i 

l,m =0 


(25) 


where 

(_ i^ k+j+l n \b a k+1 (k + 1)* T {bl + b + 1 + m) fi + j — l\fn — i 
l ’ m l\m\ (n — i)\(n — 1)! (bl + b + m)Y (bl + b + 1) \ k )\j 

and hbi+b+m is the EBS density function with power parameter bl + b + m. 
Equation (I25|) means that the density function of the WBS order statistics 
is a linear mixture of the EBS densities. Then, we can easily obtain the 
mathematical properties for T x ^ n . For example, the pth moment of T x ^ n is 

OO 

E ( : T i,n ) = {bl + l + m) T Pi(w+I+m _ 1 ). 

1 , 771 =0 


5 Estimation 


In this section, we consider estimation of the unknown parameters of the 
WBS distribution by the method of maximum likelihood. Let Xi,x 2 ,... ,x n 
be observed values of X lx X 2 ,..., X n , n independent random variables having 
the WBS distribution with unknown parameter vector £ = (a, (3, a,b) T . The 
total log-likelihood function for £, is given by 


Q lb lb 

t = t(Q=n log (a) + log (6) + log [k (a, p)\ - - ^ log (t,) + ^ log (t, + f3) 

z i =1 2=1 


i n / , \ n n 

- 2^2 S r ( ij ) + ( b “ H log M - ( b + !) X] lo § t 1 

2=1 ' ' 2=1 2=1 




$ (vj) 

.1 - ^ («i) 


$ (Vi)] 
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Then, the components of the unit score vector U = U (£) = ( d£/da , d£/d i 8, d£/da, d£/db) T 
are given by 


dt n ( _2_\ ly/ti _ (6-1) y^ Vj(t> (vi) 

da a \ a 2 ) a 3 \/3 ti) a <f> (vp 

_ (b + 1) y-> Vj(j) (Vj) ab Vj(j) (Vj) $ ( Vi) b ~ x 

M__n_ A 1 1 A/fr /A (6-l) ^ ' r (\/^)^W 

9/9 2/3 + ^ U + ,9 + 2/9a 2 ^ V/ 3 A/ 2 /9a $ (?k) 

(6 + 1 ) ™ T ( VU / t ) 0 fa ) ab ” r </ , (^) $ (^) 6 ” 1 
2 f 3a ~t 1 - $ fa) + 2^a ^ [1 - <f> (ni)] 6+1 


cM 

da 


n 

a 


£ 


_i -$fa)_ 


and 


d£ 

db 


n n 

y + Y1 log fa)l - log [1 - ^ fa)] 

i=l i =1 


-E 


$ fa) 

.1 - $fa) 


log 


$fa) 

1 - <h (vi) 


5 


where 0 (•) is the standard normal density function, r y\/ti/Pj = \JUjP + 
y/P/U and Vi = Ot 1 | \[uJP ~ fafafa} for i = 1, ..., n. The maximum 

likelihood estimate £ of £ is obtained by setting these equations to zero, 
U (£) = 0, solving them simultaneously. These equations cannot be solved 
analytically and statistical software can be used to solve them numerically 
via iterative methods such as the Newton-Raphson algorithm. 


We can use the normal approximation of the MLE of £ to construct approxi¬ 
mate confidence intervals for the parameters. Under some regular conditions 
(see [ 9 ], Chapter 9]) that are fulfilled for the parameters in the interior of the 

parameter space, the asymptotic distribution of \fn ^£ — £^ is multivariate 

normal A /4 (0, J _1 (£)), where / (£) is the expected information matrix. This 
asymptotic behavior is valid if / (£) is replaced by the observed information 

matrix, J (£), evaluated at £,i.e. J (£j . The observed information matrix is 
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given by 


J(0 


^ Laa. 

La(3 

Laa 

Lab 

\ 



Lp a 

I (3 b 




laa 

lab 


\ * 



hb 

J 


whose elements are given in the Appendix. The approximate 100(1 — rj)% 
two-sided confidence intervals for a, (3, a and b are given by a ± zr y/ var (a), 


f3 ± zr ^ jvar a ± zr y/ var (a) and b ± zr yjvar (pj respectively, where 

zr is the quantile (l — |) of the standard normal distribution and var (•) is 
the diagonal element of J _1 (M corresponding to each parameter. 


6 Applications 


In this section, we demonstrate the applicability and flexibility of the WBS 
distribution by means of two well-known real data sets with different shapes 
The first data set is given by Meeker and Escobar data [l8| and the second 


data set is given in 27]. The first data set has a bathtub shaped failure rate 


function whereas the second data set has an increasing failure rate function. 


For these data sets, we compare the results of fitting the WBS distribu¬ 
tion with the Beta BS (BBS), Kumaraswamy BS (KBS), McDonald BS 
(McBS), Marshall-Olkin extended BS (MOEBS), gamma BS (GBS), EBS 
and BS distributions using the graphical method, minus twice the maxi¬ 
mized log-likelihood (—2£), Akaike information criterion (AIC), Bayesian in¬ 
formation criterion (BIG), Consistent akaike information criterion (CAIC) 
and Kolmogorov-Smirnov (K-S) test. The pdfs of the BBS, KBS, McBS, 
MOEBS and GBS distributions (for t > 0) are given by 


h (t) = M [1-4 Mb 1 . h (t) = abg (t) 4"- 1 (v) [1 - 4“ («)]*-' , 

B (a, b) 


hit) 


eg it) 

B ( a/c , b) 


<|) a—1 


(u)[l^$ c (u)f \ 


h(x) 


_ ag it) _ 

[l-(l-a)«h(-u)] 2 


and f 5 (t) = [- log {1 - $ (u)}] 6 1 

pdf (j2J). and a, (3, a, 6, c > 0. 


respectively, where g (t) is the BS(a,/3) 
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6.1 Meeker and Escobar data 


The first data set represents failure and running times of 30 devices provided 
by Meeker and Escobar 18, p.383]. The data set is: 2, 10, 13, 23, 23, 28, 
30, 65 ,80, 88, 106, 143, 147, 173, 181, 212, 245, 247, 261, 266, 275, 293, 300, 


300, 300, 300, 300, 300, 300, 300. The total time on test (TTT) plot for the 
Meeker and Escobar data in Figure 5(a) shows a convex shape followed by 
a concave shape. This corresponds to a bathtub-shaped failure rate. Hence, 
the WBS distribution is appropriate for modeling this data set. 


Table 1 lists the MLEs and their corresponding standard errors in parentheses 
of parameters of the WBS, BBS, KBS, McBS, MOEBS, GBS, EBS and BS 
distributions for Meeker and Escobar data set. The statistics —2^, AIC, 
BIC, CAIC, K-S and its p-value are listed in Table 2 for all the distributions. 
These results show that the WBS distribution has the largest p-value and the 
smallest —2£, AIC, BIC, CAIC and K-S values. So, the WBS distribution 
gives an excellent fit than the others models for Meeker and Escobar data 
set. The histogram of this data set and the plots of the estimated densities 
of all models are given in Figure 6. From this Figure, we can conclude that 
the WBS model provides a better fit to the histogram and therefore could 
be chosen as the best model for Meeker and Escobar data. 


6.2 Turbochargers failure data 


The second data set represents the time to failure(10 3 h) of turbocharger of 


one type of engine given in Xu et al. 27]. The data set is: 1.6, 2.0, 2.6, 3.0, 


3.5, 3.9, 4.5, 4.6, 4.8, 5.0, 5.1, 5.3, 5.4, 5.6, 5.8, 6.0, 6.0, 6.1, 6.3, 6.5, 6.5, 6.7, 
7.0, 7.1, 7.3, 7.3, 7.3, 7.7, 7.7, 7.8, 7.9, 8.0, 8.1, 8.3, 8.4, 8.4, 8.5, 8.7, 8.8, 


9.0. The Figure 5(b) shows concave TTT plot for the data set, indicating 
increasing failure rate function, which can be properly accommodated by a 
WBS distribution. 


Table 3 gives the MLEs of the parameters of all models used here and their 
corresponding standard errors in parentheses. The statistics — 2£, AIC, BIC, 
CAIC, K-S and its p-value are listed in Table 4. From this Table , we 
can see the WBS distribution as the best fit for the data set among all the 
seven models. The histogram of this data set and the plots of the estimated 
densities of all models are shown in Figure 7. So, the WBS model provides 
a better fit to second data set. 


19 





Table 1: MLEs and their standard errors in parentheses for the first data. 


Model 

a 

fi 

a 

b 

c 

WBS 

0.8152 

22.9053 

0.1115 

0.2683 

— 


(0.5466) 

(13.5555) 

(0.0674) 

(0.2193) 

— 

McBS 

22.3663 

0.3293 

6.9147 

124.9055 

67.0133 


(19.2479) 

(0.5770) 

(1.3436) 

(153.8349) 

(41.2595) 

MOEBS 

1.9735 

13.8678 

17.1905 

— 

— 


(0.5232) 

(7.7663) 

(10.0909) 

— 

— 

KBS 

11.3624 

6.5795 

11.1898 

72.6776 

— 


(6.0551) 

(7.5249) 

(2.3506) 

(102.1419) 

— 

GBS 

5.6073 

1.3777 

3.6317 

— 

— 


(1.8869) 

(0.9472) 

(0.5234) 

— 

— 

BBS 

15.6640 

3.9207 

31.7249 

17.4625 

— 


(17.9647) 

(5.5777) 

(49.6707) 

(36.5502) 

— 

EBS 

4.8477 

3.8141 

5.7211 

— 

— 


(3.0942) 

(5.1989) 

(1.8078) 

— 

— 

BS 

1.6778 

64.0791 

— 

— 

— 


(0.2218) 

(14.5028) 

— 

— 

— 







Tabic 2: MLEs and their standard errors in parentheses for the second data. 


Model 

a 

P 

a 

b 

c 

WBS 

0.2007 

3.4802 

0.1185 

0.2323 

— 


(0.127081) 

(0.626594) 

(0.072865) 

(0.255233) 

— 

McBS 

10.8469 

0.0311 

21.5229 

51.3146 

59.8247 


(0.3520032) 

(0.0016003) 

(4.1772473) 

(0.1028498) 

(0.1186039) 

MOEBS 

0.5269 

2.1087 

74.3785 

— 

— 


(0.1009906) 

(0.4566437) 

(0.0049525) 

— 

— 

KBS 

7.7703 

0.1109 

23.9059 

63.4929 

— 


(2.8922) 

(0.0834) 

(3.4606) 

(0.0374) 

— 

GBS 

4.5864 

0.0160 

11.6225 

— 

— 


(0.0924) 

(0.0027) 

(1.9661) 

— 

— 

BBS 

10.9655 

0.0655 

64.5533 

15.7442 

— 


(11.8708) 

(0.1504) 

(1.9297) 

(4.1212) 

— 

EBS 

5.1071 

0.0493 

41.1747 

— 

— 


(12.3684) 

(0.2397) 

(17.0774) 

— 

— 

BS 

0.4139 

5.7538 

— 

— 

— 


(0.0463) 

(0.3684) 

— 

— 

— 



0 50 100 150 200 250 300 
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Figure 6. The histogram and the estimated densities of the first data set. 
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Table 3: The statistics: -2£, AIC, CAIC, K-S and its p-value for the first 
data._ 


Model 

—21 

AIC 

BIC 

CAIC 

K-S 

p-Value 

WBS 

352.8431 

360.8431 

366.4479 

362.4431 

0.1685 

0.3618 

McBS 

357.8659 

367.8659 

374.8719 

370.3659 

0.2273 

0.0899 

MOEBS 

363.1652 

369.1652 

373.3688 

370.0883 

0.1854 

0.2536 

KBS 

362.5005 

370.5005 

376.1053 

372.1005 

0.2054 

0.1591 

BBS 

366.5623 

374.5623 

380.1671 

376.1623 

0.2075 

0.1510 

GBS 

367.4188 

373.4188 

377.6223 

374.3418 

0.2213 

0.1058 

EBS 

368.9539 

374.9539 

379.1575 

375.8770 

0.2043 

0.1635 

BS 

385.5103 

389.5103 

392.3127 

389.9547 

0.3218 

0.0040 

Table 4: The statistics: 

-2/?, AIC, CAIC, K-S and its p-value for the second 

data. 







Model 

—21 

AIC 

BIC 

CAIC 

K-S 

p-Value 

WBS 

157.1875 

165.1875 

171.9431 

166.3304 

0.0778 

0.9685 

McBS 

164.9313 

174.9313 

183.3757 

176.696 

0.1066 

0.7535 

MOEBS 

167.0805 

173.0805 

178.1472 

173.7472 

0.0909 

0.8958 

KBS 

166.1958 

174.1958 

180.9513 

175.3387 

0.1119 

0.6976 

BBS 

173.0616 

181.0616 

187.8172 

182.2045 

0.1205 

0.6067 

GBS 

173.4768 

179.4768 

184.5434 

180.1435 

0.1199 

0.6136 

EBS 

180.8146 

186.8146 

191.8812 

187.4813 

0.1607 

0.2531 

BS 

182.7348 

186.7348 

190.1125 

187.0591 

0.1653 

0.2243 



2 4 6 8 


t 


Figure 7. The histogram and the estimated densities of the second data 

set. 
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7 Conclusions 


In this paper, we introduce a new four-parameter model called the Weibull 
Birnbaum-Saunders (WBS) distribution that extends the Birnbaum-Saunders 
distribution. The WBS hazard function can be decreasing, increasing, upside- 
down bathtub, bathtub-shaped or modified bathtub shaped depending on its 
parameters. So, the WBS model can be used quite effectively in analyzing 
lifetime data. The properties of the new distribution including expansions for 
the density function, moments, generating function, order statistics, quan¬ 
tile function, mean deviations and reliability are provided. We estimate the 
model parameters by maximum likelihood and obtain the observed informa¬ 
tion matrix. An application of the WBS distribution to two real data sets is 
used to illustrate that this distribution provides a better fit than many other 
related non-nested models. 
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8 Appendix 

The elements of the observed information matrix J (£) for the parameters 
(a, 0, a, b ) are 


j _ 6n_3- \ - / tj p \ 

aa ~ a 2 a 4 a 4 ^(/3 t I 

i= 1 v 
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